C****** NOTE DIMENSION ON ORBELM SHOULD NOT BE (1,6) BUT SHOULD
C******  BE (NREV,6)
C******  THIS IS CHANGED TO DECREASE THE CORE NEEDED FOR INVERSE
C******   RUNS. ORBELM IS USED ONLY TO DETERMINE THE EFFECT
C***** OF USING ONLY PARIAL REVS FOR USING ELEMENT RATES
C*****  FIT FOR THE GRAVITY FIELD OF THE PLANET
      PARAMETER NMAS=7001, NHARM=1, ndisk=10, NRES=721, NUMPAR=12
      PARAMETER NCDSK=10
      PARAMETER IDIM=1+5*NDISK+2*NHARM*NHARM+NHARM+4*NMAS+25*NCDSK
      PARAMETER NHDIM=NHARM+1
      PARAMETER NUM7=7*NUMPAR
      PARAMETER NREV=3
      PARAMETER NREVI=3
      PARAMETER NTRAJG=3
      PARAMETER NUMPA1=NUMPAR+1
      PARAMETER MRY=NUMPAR*NUMPA1/2+NUMPA1
C
C     NMAS = MAXIMUM NUMBER OF MASS POINTS
C     NHARM = MAXIMUM NUMBER OF HARMONICS
C     NDISK = MAXIMUM NUMBER OF DISKS
C     NRES = MAXIMUM NUMBER OF POINTS TO BE PRINTED
C     NUMPAR = MAXIMUM NUMBER OF PARAMETERS TO BE OPTIMIZED
C
      COMMON RC,XMU,COM
      COMMON /ALAT1/ LUNLAT(NMAS)
      COMMON /ALON1/ LUNLON(NMAS)
      COMMON /AMAS1/ LUNMAS(NMAS)
      COMMON /ARAD1/ LUNRAD(NMAS)
      COMMON /LABEL/ CTYP(3,10),FTYP(3,5)
      COMMON /LONG/ ELONG(NRES)
      COMMON /LAT/ ELAT(NRES)
      COMMON /ALT/ EALT(NRES)
      COMMON /ORB/ OE,T1,NOEST,NEST
      COMMON /PARA/ NOPAR,NOPTCY,TOL,NTP(NUMPAR)
      COMMON /PARB/ NTD(NUMPAR)
      COMMON /PARC/ PARDEL(NUMPAR)
      COMMON /PRAM/ KMAS,KHARM,KDISK,ICDSK
      COMMON /RESA/ ACC(NRES)
      COMMON /RESF/ FPL(NRES)
      COMMON /RESX/ NPL,XPL(NRES)
      COMMON /RESY/ YPL(NRES)
      COMMON /RESYD/ YDPL(NRES)
      COMMON /TEST/ ATABLE,LASTIT,RESSQ,DEBUG,SCPLOT,LASTPR,NCYCLE
      COMMON /TLONG/ PHLONG, PHLAT, PHALT
      COMMON /WOW/ FANOM
      COMMON /CREAM/ XSAVE(NRES)
      COMMON /PIG/ NMDATA
      COMMON /CRVDK/ AMASS,IDSWTH,NUMRN(NCDSK)
      COMMON /OPT/ NOPAR2
      CHARACTER*6  CTYP,FTYP,HELP
      INTEGER NEQ,KD(1),IFLAG,MXSTEP,KSTEP,KEMAX,KQ(6),NOEST(6),ITYPD(4)
      INTEGER PCT(4),ILAB(NUMPAR)
      REAL PHLONG,NEWDEN,ELONG, PHLAT, ELAT, PHALT
      REAL EP(6),HMINA,HMAXA,EMAX
      REAL EPZ(6)
      REAL LUNMAS,LUNRAD,LUNLAT,LUNLON
      DIMENSION W(NUMPAR), SIG(NUMPAR)
      DIMENSION XY(NRES,5), PARAM(NUMPAR), NDEX(NUMPAR)
      DIMENSION PARTAL(NRES,NUMPA1)
      DIMENSION ORBELM(1,6), ORBELP(1,6), HELP(3,15)
      DOUBLE PRECISION YDPL,XSAVE,ACC
      DIMENSION DELTA(NREVI,6), DELTP(NREVI,6)
      DOUBLE PRECISION OCCUPP
      DOUBLE PRECISION RES1(NRES),RES2(NRES),RH(3),SQ
      DOUBLE PRECISION ORBELM,ORBELP,DELTA,DELTP
      DOUBLE PRECISION DPARAM(NUMPAR)
      DOUBLE PRECISION HZ,TZ,AZ,BZ,CZ,CO,RCSQ
      DOUBLE PRECISION SAV(7,NRES),SAVFO(3,NRES)
      DOUBLE PRECISION F(6),DT(20,6),DP(6),DP1(6),J2,T1,FO(3)
      DOUBLE PRECISION T,Y(6),H,DELT,TFINAL,YN(6),EPS,EPSINV
      DOUBLE PRECISION XMU,OE(12),OEL(6,NTRAJG),DRAD,RC,R0(3),ROTA,ROTR,
     1ERS,ERS1,ERD,ROBS,OBLAT,OBLON
      DOUBLE PRECISION PNM(NHDIM,NHDIM),FF(3),CONST,COM(IDIM)
      DOUBLE PRECISION X(6),TAU,CTAU,STAU,RHO,TH,PH
      DOUBLE PRECISION AZIPP,AZIPPO
      DOUBLE PRECISION SHJ(NHDIM),SHC(NHDIM,NHDIM),SHS(NHDIM,NHDIM)
      DOUBLE PRECISION FANOM,W
      DOUBLE PRECISION PARTAL,R(MRY),DPSOS,SIG,COVR(NUMPAR,NUMPAR),WORK(
     1NUMPAR)
      DOUBLE PRECISION OBLAT1(NTRAJG),OBLON1(NTRAJG),PARAM
      LOGICAL PRINT,MNAM,INITL,PRRES,LAST,PLOTR,RESCON,PRTRAJ,TRANUR
      LOGICAL FTABLE,NUTABL,ATABLE,MORTAB,NEWCAS,NOERR,COVARI,RESPRT
      LOGICAL TOCCU,NOCCUP,NOCCU,CORLAT
      LOGICAL NOPER
      LOGICAL PRTRA1
      LOGICAL IDSWTH
      LOGICAL LASTIT,DEBUG,SCPLOT,LASTPR,DEBUG1
C      EXTERNAL PGID(FOR)
      DATA NEQ /3/, NUMPAT /NUMPAR/, KD(1) /2/, EPS /1.0D18/, EPSINV /1.
     1D-12/, CO /.1745329251994430D-1/, DRAD /57.29577951308232D0/, HELP
     2 /45*'      '/
C
      DATA (CTYP(MX,1),MX=1,3) /'     L','ATITUD','E     '/
      DATA (CTYP(MX,2),MX=1,3) /'    LO','NGITUD','E     '/
      DATA (CTYP(MX,3),MX=1,3) /'      ',' DEPTH','      '/
      DATA (CTYP(MX,4),MX=1,3) /'      ','RADIUS','      '/
      DATA (CTYP(MX,5),MX=1,3) /'  MASS',' OF RI','NG 1  '/
      DATA (CTYP(MX,6),MX=1,3) /'  MASS',' OF RI','NG 2  '/
      DATA (CTYP(MX,7),MX=1,3) /'  MASS',' OF RI','NG 3  '/
      DATA (CTYP(MX,8),MX=1,3) /'  MASS',' OF RI','NG 4  '/
      DATA (CTYP(MX,9),MX=1,3) /'  MASS',' OF RI','NG 5  '/
      DATA (CTYP(MX,10),MX=1,3) /'  MASS',' OF RI','NG 6  '/
      DATA (FTYP(MX,1),MX=1,3) /'  MASS',' OF DI','SK   '/
      DATA (FTYP(MX,2),MX=1,3) /'  RADI','US OF ','DISK  '/
      DATA (FTYP(MX,3),MX=1,3) /' LATIT','UDE OF',' DISK '/
      DATA (FTYP(MX,4),MX=1,3) /'LONGIT','UDE OF',' DISK '/
      DATA (FTYP(MX,5),MX=1,3) /' RADIA','L DIST','ANCE  '/
C
C
C  SUPPLY INITIAL CONDITIONS FOR THE INTEGRATION IN A NON-ROTATING SYSTE
C        WHOSE PLANE OF REFERENCE IS THE EQUATORIAL PLANE OF THE PRIMARY
C
C   RC     = THE EQUATORIAL RADIUS OF THE PRIMARY
C   ROTA   = THE ANGLE(MEASURED EASTWARDS IN DEGREES) THROUGH WHICH THE
C            FIXED SYSTEM MUST BE ROTATED ABOUT THE Z-AXIS TO CONCIDE
C            WITH THE SYSTEM ROTATING WITH THE PLANET AT T=0.
C   ROTR   = ROTATION RATE OF THE PLANET IN DEGREES PER SECOND
C   XMU    = GM (KM**3/SEC**2)
C   IT     = INPUT OPTION  =1  ORBITAL ELEMENTS BEING INPUT
C                          =2  POSITION AND VELOCITY INPUT
C   OEL    = EITHER A,E,M,I,W,N(ANGLES IN DEGREES)
C            OR X,Y,Z,XDOT,YDOT,ZDOT(AT TIME T)
C   T      = INITIAL TIME IN SECONDS
C   DELT   = PRINT INTERVAL
C   TFINAL = FINAL TIME
C   H      = NOMINAL STEP SIZE FOR INTEGRATION
C   HMINA  = MINIMUM STEP SIZE
C   HMAXA  = MAXIMUM STEP SIZE
C   MXSTEP = MAXIMUM NUMBER OF STEPS
C   EP     = MAXIMUM ERROR IN VELOCITY PER STEP
C   J2     = OBLATENESS FOR LONG ARC INTEGRATION
C   NCOR   = CORRECTION OPTION  =1 CORRECT FOR POSITION ONLY
C                               =2 CORRECT FOR VELOCITY ONLY
C                               =3 CORRECT FOR BOTH POSITIONS AND VELOCI
C   PRRES  = PRINT OPTION FOR DIFFERENCES AND CALCULATIONS OF THE
C            ACCELERATIONS (NOMINAL VALUE=TRUE)
C   RESPRT = PRINTS THE RESIDUALS,CALCULATED ACCELLERATIONS AND FORCES
C            (NOMINAL VALUE=.FALSE.)
C   PLOTR  = RESIDUAL PLOTS- POSITION,VELOCITY,ACCELERATION AND FORCES
C            (NOMINAL VALUE=FALSE)
C   NEST   = NUMBER OF ELEMENTS TO BE ESTIMATED
C   NOEST  = CODE TO SPECIFY WHICH ELEMENTS SHOULD BE ESTIMATED. IT IS A
C            VECTOR OF DIMENSION 6 CORRESPONDING TO THE 6 ORBITAL ELEMEN
C            FOR EACH ELEMENT  =1  ESTIMATE
C                              =0  DO NOT ESTIMATE
C   ERS    =CONVERGENCE CRITERION FOR SIGMA OF POSITION.
C   ERS1   =CONVERGENCE CRITERION FOR SIGMA OF VELOCITY.
C   ERD    =CONVERGENCE CRITERION FOR DIFFERENCE IN SIGMA FROM PREVIOUS
C           ITERATION
C   ROBS   = DISTANCE FROM THE CENTER OF THE PRIMARY TO THE OBSERVER
C   OBLAT  = LATITUDE OF THE OBSERVER
C   OBLON  = LONGITUDE OF THE OBSERVER
C   NIT    =MAXIMUM NUMBER OF ITERATIONS
C   IFIRST = INPUT OPTION
C             =1   NEW CASE READ COMPLETE DATA SET
C             =0   NO MORE DATA TO BE READ FOR THIS CASE
C   RESCON = DIVERGENT RESIDUAL CUTOFF OPTION
C             =TRUE   STOP IF RESIDUALS GET LARGER FROM ONE ITERATION TO
C                     THE NEXT
C             =FALSE   DONT TEST FOR RESIDUAL CONVERGENCE
C
C   PRTRAJ = TRAJECTORY PRINT OPTION
C   TRANUR = RUN OPTION   =TRUE - ONLY COMPUTE TRAJECTORY
C     PRINT   INPUT PRINT OPTION
C              = TRUE    PRINT DATA READ BY SUBROUTINE READS (DEFAULT
C                        VALUE)
C              = FALSE   DONT PRINT DATA READ BY SUBROUTINE READS
C     MNAM    POINT MASS READ OPTION
C              = TRUE    IF POINT MASS DATA IS TO BE READ FROM NAMELIST
C                        /MASPTS/
C              = FALSE   OTHERWISE (DEFAULT VALUE)
C     KMAS    NUMBER OF POINT MASSES (DEFAULT VALUE = NMAS)
C     KHARM   NUMBER OF SPHERICAL HARMONICS (DEFAULT VALUE = NHARM)
C     KDISK   NUMBER OF DISK MASSES (DEFAULT VALUE = NDISK)
C   FTABLE = TRUE IF FORCES ARE TO BE TAKEN FROM A TABLE, FOR POINTS
C                 WITHIN THE LONGITUDE LIMITS OF THE TABLE.
C          = FALSE OTHERWISE (DEFAULT VALUE)
C   NUTABL = TRUE ONLY IF FTABLE = TRUE AND FORCES ARE TO BE READ FROM
C                 A CARD DECK.
C          = FALSE IF FTABLE = FALSE, OR IF FTABLE = TRUE AND FORCES
C                 HAVE ALREADY BEEN READ FROM A CARD DECK ON A PREVIOUS
C                 CASE IN THIS RUN.  (DEFAULT VALUE)
C     ATABLE = TRUE IF LINE OF SIGHT ACCELERATIONS ARE TO BE TAKEN FROM
C              A TABLE FOR COMPARISONS WITH CALCULATED VALUES
C            = FALSE IF STANDARD ORBSIM RUN IS WANTED (DEFAULT VALUE)
C     MORTAB = TRUE ONLY IF ATABLE = TRUE AND A NEW SET OF ACCELERATIONS
C              IS TO BE READ FROM A CARD DECK
C            = FALSE IF ATABLE = FALSE, OR IF ATABLE = TRUE AND
C              ACCELERATIONS HAVE ALREADY BEEN READ FROM A CARD DECK ON
C              A PREVIOUS CASE IN THIS RUN (DEFAULT VALUE)
C     DEBUG  = FALSE TO AVOID MOST INTERMEDIATE PRINTOUT
C            = TRUE TO GET MOST INTERMEDIATE PRINTOUT
C     SCPLOT = TRUE IF AN SC4020 PLOT IS WANTED (NOTE - AN ASG CONTROL
C              CARD MUST BE USED TO ASSIGN THE PLOT FILE IN THIS CASE)
C            = FALSE IF NO SC4020 PLOT IS WANTED
C     NOPAR  = NUMBER OF PARAMETERS TO BE OPTIMIZED (FROM 0 THROUGH
C              NUMPAR)
C     NTP(M) = TYPE OF PARAMETER M
C              = 1 FOR MASS
C              = 2 FOR RADIUS
C              = 3 FOR LATITUDE
C              = 4 FOR LONGITUDE
C              = 5 FOR DISTANCE FROM CENTER OF PLANET
C     NTD(M) = DISK TO WHICH PARAMETER M APPLIES (MUST BE NO GREATER
C              THAN KDISK)
C     PARDEL(M) = AMOUNT BY WHICH PARAMETER M WILL BE VARIED FROM ONE
C                 ITERATION TO THE NEXT.  NOTE THAT PARDEL(1) MUST BE
C                 POSITIVE, AND THE REMAINING ONES MUST BE GREATER THAN
C                 OR EQUAL TO 0.
C     NOPTCY = NUMBER OF OPTIMIZATION CYCLES DESIRED
C     TOL    = TOLERANCE ON PARDEL(1).  WHEN THE OPTIMIZATION WANTS TO
C              CUT PARDEL(1) IN HALF, BUT FINDS THAT THE RESULT WOULD BE
C              LESS THAN THE VALUE OF TOL, THEN PHASE 1 OF THE
C              OPTIMIZATION WILL END AND PHASE 2 WILL BEGIN.
C
C     COVARI  = TRUE IF COVARIANCE MATRIX *UPPER TRIANGULAR* IS TO PRINT
C               (NOMINAL VALUE = FALSE)
C     CORLAT  = TRUE IF CORRELATION MATRIX *LOWER TRIANGULAR* IS TO PRIN
C               (NOMINAL VALUE = TRUE)
C
C INPUT   TIME IN SECONDS,   DISTANCE IN KM.,   ANGLES IN DEGREES
C
C
      NAMELIST/ELEM/RC,ROTA,ROTR,XMU,IT,OEL,T,DELT,TFINAL,H,HMINA,HMAXA
     1,MXSTEP,PE,J2,NCOR,PRRES,NOEST,NEST,ERS,ERS1,ERD,ROBS,OBLAT,OBLON
     2,NIT,IFIRST,RESCON,PRTRAJ,TRANUR,PRINT,MNAM,KMAS,KHARM,KDISK,
     3FTABLE,NUTABL,ATABLE,MORTAB,DEBUG,NOPAR1,SCPLOT,NOCCUP,NORBMX,
     4TSTART,OCCUPP,OCCULS,PRTRA1,NMTRAJ,ITRMAX,ICDSK,IDSWTH,NOPAR2,
     5COVARI,CORLAT,PLOTR,RESPRT,DEBUG1,OLDDEN,NEWDEN,EPS,EPSINV
      NAMELIST /MASPTS/ LUNMAS,LUNLAT,LUNLON,LUNRAD,ICDSK
      NAMELIST /ELEM1/ OBLAT1,OBLON1
C
C   ********************  DEBUG OUTPUT FOR PARTIAL CHECK ***************
C   ********************************************************************
C
C     BEGINNING OF PROGRAM
C      - THIS IS THE STARTING POINT FOR EACH NEW CASE READ FROM CARDS.
C      - THE PURPOSE OF THE NEXT FEW CARDS IS TO SET THE NOMINAL VALUES
C        FOR VARIOUS  INPUT ITEMS, AND TO SET THE PARAMETERS USED BY THE
C        INTEGRATION ROUTINE DVDQ.
C
C     CALL ERTRAN (17,4,PCT)
C
C      CALL PGID
      EXTERNAL HANDLER
      CALL LIB$ESTABLISH(HANDLER)
   10 NOPAR=0
      NEWDEN=1.0
      OLDDEN=1.0
      NOPER=.FALSE.
      CORLAT=.TRUE.
      COVARI=.FALSE.
      G=6.670D-8
      TOCCU=.TRUE.
      NOCCUP=.TRUE.
      PRTRAJ=.FALSE.
      TRANUR=.FALSE.
      RESCON=.TRUE.
      KMAS=NMAS
      KHARM=NHARM
      KDISK=NDISK
      PRINT=.TRUE.
      MNAM=.FALSE.
      FTABLE=.FALSE.
      NUTABL=.FALSE.
      NEWCAS=.TRUE.
      ATABLE=.FALSE.
      MORTAB=.FALSE.
      LASTIT=.FALSE.
      DEBUG=.FALSE.
      SCPLOT=.FALSE.
      ROTA=0.D0
      ROTR=0.D0
      MXSTEP=100000
      ITERTN=1
      NFILE=70
      NTRAJ=1
      NORB=0
      IFIRST=1
      PRTRA1=.TRUE.
      IPARTL=0
      KGO=1
      FILMOD=0.0
      NOPAR1=0
      NOPAR2=0
      PRRES=.TRUE.
      PLOTR=.FALSE.
      RESPRT= .FALSE.
C
C     READ NEW DATA
C      - THE NEXT SECTION READS NAMELIST ELEM, AND DEPENDING ON VALUES
C        FOUND THERE, READS OTHER INPUT DATA.
C      - MANY OF THE INPUT VALUES ARE CHECKED FOR VALIDITY.  IF INVALID
C        DATA IS FOUND, THE CASE IS SKIPPED AND THE NEXT CASE IS SOUGHT.
C
      READ (5,ELEM,END=1070)
      EP(1)=PE
      NOPAR=NOPAR1+NOPAR2
C
      DO 20 IRL=1,NOPAR
        ILAB(IRL)=IRL
   20   CONTINUE
C
      AMASS=XMU*1.0D15/G
      OBLAT1(1)=OBLAT
      OBLON1(1)=OBLON
      IF (NMTRAJ.GT.1) GO TO 1060
   30 WRITE (6,1080)
      NOCCU=NOCCUP
      IF (.NOT.NOCCU) TOCCU=.FALSE.
      WRITE (6,ELEM)
      IF(NMTRAJ.GT.1) WRITE(6,ELEM1)
      DELSAV=DELT
      DELSA1=DELT
      RCSQ=RC**2
      LASTPR=LASTIT.OR.DEBUG
      IF (.NOT.FTABLE) NUTABL=.FALSE.
      IF (.NOT.ATABLE) MORTAB=.FALSE.
      IF (.NOT.NEWCAS) GO TO 40
      NEWCAS=.FALSE.
      IF (IFIRST.EQ.0) GO TO 40
      IFIRST=0
      IF (MNAM) READ (5,MASPTS)
      CALL READS (PRINT,MNAM,SHJ,SHC,SHS,NHDIM)
      IF (ICDSK.GT.0) CALL RDSK1 (NOPER,PRINT)
      IF (NOPER) GO TO 10
      IF (MORTAB) CALL INADCK (*10)
      IF (NUTABL) CALL INFDCK (*10,NEWDEN,OLDDEN)
   40 IF (NOPAR.LT.0) NOPAR=0
      IF (NOPAR.LE.NUMPAR) GO TO 50
      WRITE (6,1090) NOPAR,NUMPAT
      GO TO 10
C
C     STORE INITIAL VALUES SO THAT THE PROGRAM CAN BE REINITIALIZED FOR
C     SUCCEEDING ITERATIONS DURING THE OPTIMIZATION.
C
   50 DO 60 I=1,6
        EPZ(I)=EP(I)
   60   CONTINUE
      DIFFER=0
      DIFFEL=0
      HZ=H
      HMINAZ=HMINA
      NCORZ=NCOR
      NESTZ=NEST
      TZ=T
      KPRINT=1
      N1=NMAS
      N2=NHARM
      N3=NDISK
      IF (KMAS.GT.N1) GO TO 70
      IF (KHARM.GT.N2) GO TO 70
      IF (KDISK.LE.N3) GO TO 80
   70 WRITE (6,1100) KMAS,KHARM,KDISK,N1,N2,N3
      GO TO 10
C
C ITYPD(I) TELLS WHICH KINDS OF DATA WERE READ -
C      (1) IS NONZERO IF MASPTS,    0 IF NOT
C      (2) IS NONZERO IF HARMONICS, 0 IF NOT
C      (3) IS NONZERO IF DISKS,     0 IF NOT
C
   80 ITYPD(1)=KMAS
      ITYPD(2)=KHARM
      ITYPD(3)=KDISK
      ITYPD(4)=ICDSK
C      MAXTIM=PCT(4)-(3600*KDISK+32700)*KDISK
      NTEST=(TFINAL-T)/DELT+0.999D0
      IF (NTEST.LT.NRES) GO TO 90
      NRESA=NRES
      WRITE (6,1110) T,TFINAL,DELT,NTEST,NRESA
      GO TO 10
   90 CONTINUE
C
C     OPTIMIZATION CONTROL
C      - THIS SECTION CHECKS THE OPTIMIZATION PARAMETERS, IF VARIABLE
C        NOPAR IS POSITIVE.
C      - OTHERWISE, IT SETS UP TO PERFORM JUST ONE CALCULATION FOR THE
C        CASE.
C
      ASSIGN 200 TO NEXT
      NCYCLE=0
      IF (NOPAR.EQ.0) GO TO 410
      ASSIGN 170 TO NEXT
      NOERR=.TRUE.
      DO 100 M=1,NOPAR
        IF (NTP(M).LT.1) NOERR=.FALSE.
        IF (NTD(M).LT.1) NOERR=.FALSE.
        IZZZZZ=MAX0(KDISK,ICDSK)
        IF (NTD(M).GT.IZZZZZ) NOERR=.FALSE.
  100   CONTINUE
      IF (NOERR) GO TO 110
      WRITE (6,1120) KDISK,(M,NTP(M),NTD(M),PARDEL(M),M=1,NOPAR)
      GO TO 10
  110 ID=4*KMAS+2*KHARM**2+KHARM-KDISK
      IF (NOPAR1.EQ.0) GO TO 130
      DO 120 M=1,NOPAR
        IM=ID+NTD(M)+KDISK*NTP(M)
        NDEX(M)=IM
        DPARAM(M)=COM(IM)
        PARAM(M)=0.0
  120   CONTINUE
  130 IF (NOPAR2.EQ.0) GO TO 370
      ID=ID+6*KDISK
      NOPAR3=NOPAR1+1
      DO 160 M=NOPAR3,NOPAR
        ITRASH=NTD(M)-1
CHANGE MADE: 10/02/84 ... R.N. WIMBERLY  -INFORMATICS GENERAL CORPORATION-
C... CHANGE- IM=0  TO  IM=ID  IN ORDER TO OBTAIN THE APPROPRIATE STARTING
C    LOCATION IN THE COM ARRAY
C
        IM=ID
        IF (ITRASH.EQ.0) GO TO 150
        DO 140 MM=1,ITRASH
          IM=4+NUMRN(MM)+IM
  140     CONTINUE
  150   IM=IM+NTP(M)
C
CHANGE MADE: 10/02/84 ...
C...  CHANGED- NDEX(M)=IM+ID  TO  NDEX(M)=IM  FOR OBVIOUS REASONS
C
        NDEX(M)=IM
        DPARAM(M)=COM(IM)
        PARAM(M)=0.0
  160   CONTINUE
      GO TO 370
C
C     AFTER EACH STEP IN THE OPTIMIZATION (WHICH INVOLVES A COMPLETE
C     TRAJECTORY CALCULATION PLUS A FEW ITERATIONS), COME TO THIS POINT
C     TO CALCULATE THE RESIDUAL SUM OF SQUARES, FOR USE IN THE DIRECT
C     SEARCH ROUTINE.
C
  170 IF (KGO.GT.1) GO TO 200
      IF (IPARTL.GT.NOPAR) GO TO 180
      IM=NDEX(IPARTL)
      PARAM(IPARTL)=PARDEL(IPARTL)
      GO TO 420
  180 DO 190 I=1,NMDATA
  190   PARTAL(I,NOPAR+1)=-RES1(I)
      NCS1=NOPAR+1
      NOSTOC=0
C   *************************************************************
      CALL THH (PARTAL,NRES,NCS1,R,MRY,NOPAR,NOSTOC,NMDATA,DPSOS,FILMOD)
C
  200 IPARTL=0
      NTRAJ=NTRAJ+1
      IF (NTRAJ.GT.NMTRAJ) GO TO 210
      OBLAT=OBLAT1(NTRAJ)
      OBLON=OBLON1(NTRAJ)
      IF (NOPAR.EQ.0) GO TO 440
      CALL INADCK (*10)
      CALL INFDCK(*10,NEWDEN,OLDDEN)
      GO TO 420
  210 IF (NOPAR.EQ.0) GO TO 10
      IF (KGO.GT.1) GO TO 320
      CALL UPINV (R,NOPAR,R,W)
      SIG1=0.
      NRM=NOPAR*(NOPAR+1)/2
      JJ=0
      DO 230 I=1,NOPAR
        JJ=JJ+I
        J=JJ
        PARAM(I)=0.D0
        SIG(I)=0.D0
        DO 220 KK=I,NOPAR
          PARAM(I)=PARAM(I)+R(J)*R(NRM+KK)
          SIG(I)=SIG(I)+R(J)**2
  220     J=J+KK
  230   SIG(I)=DSQRT(SIG(I))
C   ****************  DEBUG  PRINT TO CHECK PARTIALS  ***********
C   *************************************************************
      DO 240 M=1,NOPAR
        IM=NDEX(M)
        DPARAM(M)=DPARAM(M)+PARAM(M)
        COM(IM)=DPARAM(M)
  240   PARAM(M)=0.0
      ITERTN=ITERTN+1
      KGO=1
      IF (ITERTN.GT.NOPTCY) KGO=2
      IF (RESSQ1.LT.1.0D-8) KGO=3
      IF (KGO.LT.3) RESSQ1=0
      NTRAJ=1
      IF (MORTAB) CALL INADCK (*10)
      OBLAT=OBLAT1(1)
      OBLON=OBLON1(1)
      WRITE (6,1130)
C
      IF (NOPAR1.LT.1) GO TO 270
      DO 260 ICHK=1,NOPAR1
        INDX=NTP(ICHK)
        DO 250 MHP=1,3
          HELP(MHP,ICHK)=FTYP(MHP,IDNX)
  250     CONTINUE
  260   CONTINUE
  270 CONTINUE
      IF (NOPAR2.LT.1) GO TO 300
      DO 290 ICHK=1,NOPAR2
        ICK=NOPAR1+ICHK
        IDNX=NTP(ICK)
        DO 280 MHP=1,3
          HELP(MHP,ICK)=CTYP(MHP,IDNX)
  280     CONTINUE
  290   CONTINUE
  300 CONTINUE
C
      DO 310 M=1,NOPAR
        IM=NDEX(M)
        WRITE (6,1140) (HELP(MX,M),MX=1,3),M,COM(IM),SIG(M)
  310   CONTINUE
C
CALL SUBROUTINE CORREL TO COMPUTE THE CORRELATION MATRIX (LOWER TRIANGUL
C     ******************************************************************
      IF (CORLAT) CALL CORREL (R,NOPAR,SIG,ILAB,WORK)
CALL SUBROUTINE COVAR TO COMPUTE THE COVARIANCE MATRIX (UPPER TRIANGULAR
      IF (COVARI) CALL COVAR (NOPAR,R,ILAB,COVR)
C     ******************************************************************
      IF (.NOT.LASTIT) GO TO 370
  320 IF (KDISK.EQ.0) GO TO 340
      WRITE (6,1150)
      DO 330 I=1,KDISK
        I1=ID+I+KDISK
        I2=I1+KDISK
        I3=I2+KDISK
        I4=I3+KDISK
        I5=I4+KDISK
        WRITE (6,1160) I,COM(I1),COM(I2),COM(I3),COM(I4),COM(I5)
  330   CONTINUE
  340 ID=ID+1
      DO 360 I=1,ICDSK
        ILAT=ID
        ILON=ID+1
        IRADC=ID+2
        IRNGD=ID+3
        WRITE (6,1170) I,COM(ILAT),COM(ILON),COM(IRADC),COM(IRNGD)
        IRNUM=NUMRN(I)
        DO 350 IBLIB=1,IRNUM
          ICRAM=IRNGD+IBLIB
          WRITE (6,1180) IBLIB,COM(ICRAM)
  350     CONTINUE
        ID=1+ICRAM
  360   CONTINUE
      GO TO 10
  370 CONTINUE
      GO TO (420,380,390), KGO
  380 WRITE (6,1190)
      GO TO 410
  390 WRITE (6,1200)
      GO TO 410
  400 CONTINUE
  410 LASTIT=.TRUE.
      LASTPR=.TRUE.
      IF (NOPAR.EQ.0) GO TO 440
C
C     CALCULATION OF A SINGLE CASE
C      - THIS SECTION FIRST REINITIALIZES CERTAIN VARIABLES.
C      - THEN IT INITIALIZES THE INTEGRATION ROUTINE DVDQ.
C      - THEN IT INTEGRATES THE BASIC DIFFERENTIAL EQUATIONS.
C      - THEN IT ITERATES ON CERTAIN PARAMETERS (WITHIN A SINGLE
C        OPTIMIZATION STEP).
C
  420 DO 430 M=1,NOPAR
        IM=NDEX(M)
        COM(IM)=DPARAM(M)+DBLE(PARAM(M))
  430   CONTINUE
  440 T1=TZ
      DO 450 I=1,6
        EP(I)=EPZ(I)
  450   CONTINUE
      AZ=OBLAT*CO
      BZ=OBLON*CO
      CZ=ROBS*DCOS(AZ)
      R0(1)=CZ*DCOS(BZ)
      R0(2)=CZ*DSIN(BZ)
      R0(3)=ROBS*DSIN(AZ)
      NCYCLE=NCYCLE+1
      H=HZ
      HMINA=HMINAZ
      NCOR=NCORZ
      NEST=NESTZ
      T=TZ
      GO TO (460,470), IT
C
C  POS AND VEL FROM ORBITAL ELEMENTS
C  OEL = A,E,M,I,W,N    ALL ANGLES IN DEGREES.
C
  460 OE(1)=OEL(1,NTRAJ)
      OE(2)=OEL(2,NTRAJ)
      OE(3)=OEL(3,NTRAJ)*CO
      AZ=OEL(4,NTRAJ)*CO
      OE(4)=DSIN(AZ)
      OE(5)=DCOS(AZ)
      AZ=OEL(5,NTRAJ)*CO
      OE(6)=DSIN(AZ)
      OE(7)=DCOS(AZ)
      AZ=OEL(6,NTRAJ)*CO
      OE(8)=DSIN(AZ)
      OE(9)=DCOS(AZ)
      OE(10)=0.D0
      CALL OSCUL (1,T,XMU,Y,OE)
      GO TO 490
C
C  OEL = X,Y,Z,XDOT,YDOT,ZDOT
C
  470 DO 480 I=1,3
        Y(I+I-1)=OEL(I,NTRAJ)
  480   Y(I+I)=OEL(I+3,NTRAJ)
      CALL OSCUL (2,T,XMU,Y,OE)
  490 CONTINUE
C
C
C  INITIALIZE DVDQ
C
      PERIOD=OE(1)*2.*3.14159*SQRT(OE(1)/XMU)
      PERIOD=7200.
      CALL DVDQ (NEQ,T,Y,F,KD,EP,IFLAG,H,HMINA,HMAXA,DELT,TFINAL,MXSTEP,
     1KSTEP,KEMAX,EMAX,KQ,YN,DT)
      INITL=.TRUE.
      N=0
      DO 500 I=1,6
        DP(I)=0.D0
  500   DP1(I)=0.D0
      GO TO 520
C
C  ENTER HERE FOR A NEW INTEGRATION STEP
C
  510 CALL DVDQ1(NEQ,KD,IFLAG,MXSTEP,KSTEP,KEMAX,KQ,HMINA,HMAXA,
     1           EMAX,T,Y,F,H,DELT,TFINAL,YN,DT,EP)
  520 GO TO (530,530,610,680,680,1030,1040,1050), IFLAG
C
C  COMPUTE THE FORCES
C
C
CONVERT TO THE ROTATING SYSTEM FIXED IN THE PLANET
C
  530 TAU=(ROTA+ROTR*T)*CO
      STAU=DSIN(TAU)
      CTAU=DCOS(TAU)
      X(1)=Y(1)*CTAU+Y(3)*STAU
      X(2)=-Y(1)*STAU+Y(3)*CTAU
      X(3)=Y(5)
      RHO=DSQRT(Y(1)*Y(1)+Y(3)*Y(3)+Y(5)*Y(5))
      TH=DASIN(X(3)/RHO)*DRAD
      PH=DATAN2(X(2),X(1))*DRAD
      FF(1)=0.0D0
      FF(2)=0.0D0
      FF(3)=0.0D0
      IF (FTABLE) CALL FCALC (*10,PH,FF)
      DO 600 I=1,4
        IF (ITYPD(I).NE.0) GO TO (540,550,570,560), I
        GO TO 600
540   IF(KMAS.EQ.3500 .AND. (RHO-RC).GT.10000.0D0) GO TO 600
        CALL FUNCT (RHO,TH,PH,EPS,FO)
        GO TO 580
  550   CALL LSHARM (RHO,TH,PH,FO,NHDIM,PNM,SHJ,SHC,SHS)
        GO TO 580
  560   CALL FCDSK (RHO,TH,PH,EPS,FO)
        GO TO 580
  570   CALL DISKS (RHO,TH,PH,EPS,FO)
  580   DO 590 IX=1,3
          FF(IX)=FO(IX)+FF(IX)
  590     CONTINUE
  600   CONTINUE
C
C  COMPUTE FORCES IN INERTIAL SYSTEM
C
      CONST=-XMU/RHO**3
      FO(1)=FF(1)*CTAU-FF(2)*STAU
      FO(2)=FF(1)*STAU+FF(2)*CTAU
      FO(3)=FF(3)
      F(1)=FO(1)+CONST*Y(1)
      F(2)=FO(2)+CONST*Y(3)
      F(3)=FO(3)+CONST*Y(5)
      GO TO 510
C
C  OUTPUT UNDER CONTROL OF DELT
C
  610 IF (NOCCUP) GO TO 630
      DELT=DELSA1
      IF (T.LT.TSTART.OR.TOCCU) GO TO 620
      IF (NPL.GT.NRES) WRITE (6,1220)
      GO TO 680
  620 DIFFEL=DIFFER
      TTEST=T+DELT
      IF (TOCCU.OR.TTEST.LE.TSTART) GO TO 630
      DELSA1=TTEST-TSTART
      DELT=TSTART-T
  630 RHO=Y(1)**2+Y(3)**2+Y(5)**2
      CONST=-XMU/(RHO*DSQRT(RHO))
      FO(1)=F(1)-CONST*Y(1)
      FO(2)=F(2)-CONST*Y(3)
      FO(3)=F(3)-CONST*Y(5)
      IF (DABS(FO(1))+DABS(FO(2))+DABS(FO(3)).GT.EPSINV) GO TO 640
      FO(1)=0.0D0
      FO(2)=0.0D0
      FO(3)=0.0D0
  640 IF (IPARTL.NE.0) GO TO 650
      CALL OUTP (ROTA,ROTR,T,Y,RC,XMU,PRTRA1)
      CALL CORECT (XMU,Y,J2,RC,T,R0,DP,NCOR,INITL,.FALSE.,PRTRA1,FO)
C
C  NCOR= 1  CORRECT POSITION ONLY
C        2  CORRECT VELOCITY ONLY
C        3  CORRECT BOTH
C
  650 INITL=.FALSE.
      N=N+1
      IF (NOCCUP) GO TO 660
      DIFFER=(FANOM-OCCUPP)
      IF (ABS(DIFFER).GT.5.) GO TO 660
      IF (DIFFER.LE.0) GO TO 660
      IF (DIFFEL.GT.0) GO TO 660
      INPNT=N
      IF (ABS(DIFFEL).LT.ABS(DIFFER)) INPNT=N-1
  660 SAV(1,N)=T
      DO 670 I=1,6
  670   SAV(I+1,N)=Y(I)
      SAVFO(1,N)=FO(1)
      SAVFO(2,N)=FO(2)
      SAVFO(3,N)=FO(3)
      GO TO 510
C
C  T=TFINAL(IFLAG=4) OR KSTEP=MXSTEP(IFLAG=5)
C
  680 RHO=Y(1)**2+Y(3)**2+Y(5)**2
      CONST=-XMU/(RHO*DSQRT(RHO))
      FO(1)=F(1)-CONST*Y(1)
      DELT=DELSA1
      FO(2)=F(2)-CONST*Y(3)
      FO(3)=F(3)-CONST*Y(5)
      IF (DABS(FO(1))+DABS(FO(2))+DABS(FO(3)).GT.EPSINV) GO TO 690
      FO(1)=0.0D0
      FO(2)=0.0D0
      FO(3)=0.0D0
  690 IF (IPARTL.NE.0) GO TO 700
      CALL OUTP (ROTA,ROTR,T,Y,RC,XMU,PRTRA1)
      CALL CORECT (XMU,Y,J2,RC,T,R0,DP,NCOR,.FALSE.,.TRUE.,PRTRA1,FO)
  700 N=N+1
      SAV(1,N)=T
      DO 710 I=1,6
        SAV(I+1,N)=Y(I)
  710   CONTINUE
      SAVFO(1,N)=FO(1)
      SAVFO(2,N)=FO(2)
      SAVFO(3,N)=FO(3)
      IF (LASTPR.AND.NOCCU) WRITE (6,1230) KSTEP,EP(1),H,IFLAG
      IF (NPL.EQ.0.OR.TRANUR) GO TO NEXT, (10,170)
      CALL PRINTR (NCOR,.FALSE.,ERS,ERS1,ERD,.FALSE.,XY,NRES,
     1 .FALSE.,.FALSE.)
      IF (IPARTL.EQ.0) GO TO 770
C***********************************C
      DO 727 ILOOP=1,2
      IF(IPARTL.GT.0) LASTPR= .FALSE.
      LAST=.FALSE.
      INITL=.TRUE.
      DO 723 I=1,6
         DP(I)=DP(I)+DP1(I)
         DP1(I)=DP(I)
  723 CONTINUE
      DO 725 I=1,N
         T=SAV(1,I)
      DO 724 J=1,6
  724    Y(J)=SAV(J+1,I)
      FO(1)=SAVFO(1,I)
      FO(2)=SAVFO(2,I)
      FO(3)=SAVFO(3,I)
      CALL OUTP(ROTA,ROTR,T,Y,RC,XMU,.FALSE.)
        IF(I.EQ.N) LAST=.TRUE.
      CALL CORECT(XMU,Y,J2,RC,T,R0,DP,NCOR,INITL,LAST,.FALSE.,FO)
      INITL=.FALSE.
C ********************** C
  725 CONTINUE
C
  726 CONTINUE
  727 CONTINUE
C
  731 CONTINUE
      OE(3)=OE(3)*DRAD
      OE(4)=DATAN2(OE(4),OE(5))*DRAD
      OE(5)=DATAN2(OE(6),OE(7))*DRAD
      OE(6)=DATAN2(OE(8),OE(9))*DRAD
      DO 732  MX=1,6
      DP(MX)=DP(MX)+DP1(MX)
      IF(MX.GT.2) DP(MX)=DP(MX)*DRAD
      OE(MX)=OE(MX)+DP(MX)
  732 CONTINUE
C*********************************C
      DO 740 I=1,N
        YDPL(I)=0.0
      SQ=0.0
      DO 720 J=1,3
      RH(J)=R0(J)-SAV(2*J,I)
  720 SQ=SQ+RH(J)*RH(J)
      SQ=DSQRT(SQ)
      DO 730 J=1,3
      RH(J)=RH(J)/SQ
  730 YDPL(I)=RH(J)*SAV(2*J+1,I)+YDPL(I)
  740 YDPL(I)=(YDPL(I)-XSAVE(I))*1.0D5
C ***************************MORE DEBUG OUTPUT*****************************
      IF(.NOT.DEBUG1) GO TO 748
      WRITE(6,741)
  741 FORMAT(///,'   ***   PARTIAL DEBUG PRINT   ***')
      DO 742 MMMM=1,3
      IF(MMMM.EQ.1) WRITE(6,743)
      IF(MMMM.EQ.2) WRITE(6,744)
      IF(MMMM.EQ.3) WRITE(6,745)
      WRITE(6,746) (SAV(2*MMMM,NR),NR=1,N)
      WRITE(6,747) (SAV(2*MMMM+1,NR),NR=1,N)
  742 CONTINUE
  743 FORMAT(1H ,'** X **')
  744 FORMAT(1H ,'** Y **')
  745 FORMAT(1H ,'** Z **')
  746 FORMAT(1H ,'  POSITION - ',5(D15.8,1X))
  747 FORMAT(1H ,'  VELOCITY - ',5(D15.8,1X))
  748 CONTINUE
C *************************************************************************
      DX=(XPL(2)-XPL(1))*6.
      NZIPP=N-1
      DO 750 I=2,NZIPP
        ACC(I)=(YDPL(I+1)-YDPL(I-1))/(2.*DX)
  750   CONTINUE
      ACC(1)=(-3.0*YDPL(1)+4.0*YDPL(2)-YDPL(3))/(2.*DX)
      ACC(N)=(3.00*YDPL(N)-4.0*YDPL(N-1)+YDPL(N-2))/(2.*DX)
C *************************************************************************
C *************************************************************************
      CALL ACALC (RES2,.TRUE.)
      DO 760 I=1,NMDATA
C *************************************************************************
      IF(DEBUG1) WRITE(6,1400)  RES2(I),RES1(I)
 1400 FORMAT(1H ,' *****  RES2=',D15.8,'   RES1=',D15.8,' *****')
C *************************************************************************
  760   PARTAL(I,IPARTL)=(RES2(I)-RES1(I))/PARAM(IPARTL)
      PARAM(IPARTL)=0.0
      IPARTL=IPARTL+1
      GO TO 170
  770 DO 850 ITER=1,NIT
        DO 780 I=1,6
          IF (DABS(DP(I)).GT.1.D-8) GO TO 790
  780     CONTINUE
        GO TO 830
  790   LAST=.FALSE.
        INITL=.TRUE.
        DO 800 I=1,6
          DP(I)=DP(I)+DP1(I)
  800     DP1(I)=DP(I)
        DO 820 I=1,N
          T=SAV(1,I)
          DO 810 J=1,6
  810       Y(J)=SAV(J+1,I)
          FO(1)=SAVFO(1,I)
          FO(2)=SAVFO(2,I)
          FO(3)=SAVFO(3,I)
          CALL OUTP (ROTA,ROTR,T,Y,RC,XMU,PRTRAJ)
          IF (I.EQ.N) LAST=.TRUE.
          CALL CORECT (XMU,Y,J2,RC,T,R0,DP,NCOR,INITL,LAST,PRTRAJ,FO)
          INITL=.FALSE.
  820     CONTINUE
        IF (NPL.EQ.0) GO TO NEXT, (10,170)
        IF (LASTPR.AND.NOCCU) WRITE (6,1240) DP
        IF (ITER.NE.NIT) GO TO 850
  830   CONTINUE
  840  CALL PRINTR(NCOR,PLOTR,ERS,ERS1,ERD,RESCON,XY,
     1 NRES,PRRES,RESPRT)
       IF(PLOTR) GO TO 860
  850   CONTINUE
  860 OE(3)=OE(3)*DRAD
      OE(4)=DATAN2(OE(4),OE(5))*DRAD
      OE(5)=DATAN2(OE(6),OE(7))*DRAD
      OE(6)=DATAN2(OE(8),OE(9))*DRAD
      IF (LASTPR.AND.NOCCU) WRITE (6,1250) (OE(I),I=1,6)
      DO 870 I=1,6
        DP(I)=DP(I)+DP1(I)
        IF (I.GT.2) DP(I)=DP(I)*DRAD
        OE(I)=OE(I)+DP(I)
  870   CONTINUE
      IF (.NOT.NOCCUP) GO TO 930
      IF (LASTPR) WRITE (6,1260) (OE(I),I=1,6)
      IF((NOPAR .EQ. 0) .AND. .NOT.(ATABLE)) GO TO 200
      CALL ACALC (RES1,.FALSE.)
      IF (NOPAR .EQ. 0) GO TO 200
      IPARTL=IPARTL+1
      RESSQ1=RESSQ1+RESSQ
      GO TO NEXT, (200,170)
  880 TFINAL=TZ+PERIOD
      NORB=NORB+1
      DO 890 I=1,6
        ORBELP(NORB,I)=OE(I)
        OE(I)=OE(I)-DP(I)
        DP1(I)=0
  890   DP(I)=0
      OE(3)=OE(3)*CO
      AZ=OE(4)*CO
      AZIPP=OE(5)*CO
      AZIPPO=OE(6)*CO
      OE(4)=DSIN(AZ)
      OE(5)=DCOS(AZ)
      OE(6)=DSIN(AZIPP)
      OE(7)=DCOS(AZIPP)
      OE(8)=DSIN(AZIPPO)
      OE(9)=DCOS(AZIPPO)
      INITL=.TRUE.
      LAST=.FALSE.
      N=N-1
      DO 910 I=1,N
        T=SAV(1,I)
        DO 900 J=1,6
  900     Y(J)=SAV(J+1,I)
        FO(1)=SAVFO(1,I)
        FO(2)=SAVFO(2,I)
        FO(3)=SAVFO(3,I)
        CALL OUTP (ROTA,ROTR,T,Y,RC,XMU,PRTRAJ)
        CALL CORECT (XMU,Y,J2,RC,T,R0,DP,NCOR,INITL,LAST,PRTRAJ,FO)
        INITL=.FALSE.
        TOCCU=.TRUE.
  910   CONTINUE
      IZIPPO=N+1
      T=SAV(1,IZIPPO)
      DO 920 J=1,6
  920   Y(J)=SAV(J+1,IZIPPO)
      DELSA1=DELSAV
      GO TO 510
  930 IF (.NOT.TOCCU) GO TO 880
      WRITE (6,1240) DP
      WRITE (6,1160) N
      DO 940 I=1,6
        ORBELM(NORB,I)=OE(I)
  940   EP(I)=EPZ(I)
      DO 950 I=1,3
        OEL(I,NTRAJ)=SAV(I+I,INPNT)
        OEL(I+3,NTRAJ)=SAV(I+I+1,INPNT)
  950   CONTINUE
      IT=2
      TOCCU=.FALSE.
      IF (NORB.GE.NORBMX) GO TO 970
      DO 960 I=1,6
        EP(I)=EPZ(I)
  960   CONTINUE
      HMINA=HMINAZ
      NCOR=NCORZ
      NEST=NESTZ
      T=SAV(1,INPNT)
      TSTART=T+OCCULS
      H=HZ
      TFINAL=T+PERIOD
      GO TO 50
  970 WRITE (6,1270)
      WRITE (6,1280)
      DO 980 NORBI=1,NORBMX
        WRITE (6,1290) (ORBELM(NORBI,I),I=1,6)
  980   CONTINUE
      WRITE (6,1300)
      WRITE (6,1280)
      DO 1000 NORBI=1,NORBMX
        WRITE (6,1290) (ORBELP(NORBI,I),I=1,6)
        IF (NORBI.EQ.NORBMX) GO TO 1000
        DO 990 I=1,6
          DELTA(NORBI,I)=ORBELM(NORBI+1,I)-ORBELM(NORBI,I)
          DELTP(NORBI,I)=ORBELP(NORBI+1,I)-ORBELP(NORBI,I)
  990     CONTINUE
 1000   CONTINUE
      NORBMX=NORBMX-1
      WRITE (6,1310)
      WRITE (6,1270)
      DO 1010 NORBI=1,NORBMX
        WRITE (6,1290) (DELTA(NORBI,I),I=1,6)
 1010   CONTINUE
      WRITE (6,1300)
      DO 1020 NORBI=1,NORBMX
        WRITE (6,1290) (DELTP(NORBI,I),I=1,6)
 1020   CONTINUE
      NORBMX=NORBMX+1
      GO TO 10
C
C  INCREASE EP TO REDUCE ROUND-OFF
C
 1030 EP(1)=32.*EMAX*EP(1)
      WRITE (6,1320) EP(1),T
      GO TO 510
C
C  ABS(H) NEEDS TO BE LESS THAN HMINA TO MEET REQUIRED ACCURACY
C
 1040 HMINA=DABS(H)
      WRITE (6,1330) H
      GO TO 510
 1050 WRITE (6,1340)
      GO TO NEXT, (10,170)
 1060 READ (5,ELEM1)
      OBLAT=OBLAT1(1)
      OBLON=OBLON1(1)
      GO TO 30
 1070 STOP
C
 1080 FORMAT (1H1)
 1090 FORMAT (9H1NOPAR = ,I4,35H, EXCEEDING MAXIMUM VALUE NUMPAR = ,I4)
 1100 FORMAT (22H1KMAS, KHARM, KDISK = ,3I4,20H, EXCEEDING MAXIMUM ,28HV
     1ALUES NMAS, NHARM, NDISK = ,3I4)
 1110 FORMAT (30H1INPUT DATA GIVES INITIAL T = ,D12.6,12H, FINAL T = ,D1
     12.6,23H, AND PRINT INTERVAL = ,D12.6/22H THIS REQUIRES SAVING ,I3,
     222H POINTS, WHICH EXCEEDS,16H THE MAXIMUM OF ,I3,26H ALLOWED BY PA
     3RAMETER NRES)
 1120 FORMAT (36H1OPTIMIZATION PARAMETER OUT OF RANGE/10H PARAMETER,5X,1
     1HM,5X,3HNTP,5X,3HNTD,10X,6HPARDEL/6H RANGE,15X,3H1-5,4X,2H1-,I2,10
     2X,27H.GE. 0   (PARDEL(1) .GT. 0)/(I16,2I8,E16.8))
 1130 FORMAT (1H0,'PARAMETER VALUES FOR ABOVE RESULTS',//,7X,'TYPE',8X,'
     1NUMBER',4X,'PARAMETER',12X,'SIGMA')
 1140 FORMAT (1X,3A6,I4,3X,D15.8,3X,D15.8)
 1150 FORMAT (16H1FINAL DISK DATA//6H  DISK,6X,9HDISK MASS,15H    DISK R
     1ADIUS,7X,8HLATITUDE,6X,20HLONGITUDE  DISTANCE ,4HFROM/70X,11HMOON
     2CENTER/)
 1160 FORMAT (I6,5D15.8)
 1170 FORMAT (1X,I5,13X,7D15.8)
 1180 FORMAT (20X,I5,6X,D12.6)
 1190 FORMAT (30H1CONVERGENCE TERMINATED AFTER ,27HMAXIMUM ALLOWED ITERA
     1TIONS.)
 1200 FORMAT (30H1CONVERGENCE TERMINATED AFTER ,
     126HSATISFYING TOLERANCE TEST.)
 1210 FORMAT (46H1CONVERGENCE TERMINATED BY MAXIMUM TIME LIMIT.)
 1220 FORMAT (1X,'NPL IS GREATER THAN NRES TURKEY')
 1230 FORMAT (1H0,'STEPS=',I5,'    ERROR=',E10.3,'    H=',D12.4,'   IFLA
     1G=',I3)
 1240 FORMAT (4H0DP=6D20.8)
 1250 FORMAT (1H1'INPUT ORBITAL ELEMENTS'/3H A=F10.4,2X'E='F6.4,2X'M='F1
     11.6,2X'I='F11.6,2X'W='F11.6,2X'N='F11.6//)
 1260 FORMAT (1H0'CORRECTED ORBITAL ELEMENTS'/3H A=F10.4,2X'E='F6.4,2X'M
     1='F11.6,2X'I='F11.6,2X'W='F11.6,2X'N='F11.6//)
 1270 FORMAT (1H1,'FULL ORBIT INTEGRATION')
 1280 FORMAT (1X,4X,'A',7X,'E',12X,'M',13X,'I',13X,'W',13X,'N')
 1290 FORMAT (1X,F10.4,2X,F8.6,2X,F11.6,2X,F11.6,2X,F11.6,2X,F11.6)
 1300 FORMAT (1H1,'PARTIAL ORBIT INTEGRATION')
 1310 FORMAT (1X,'DELTAS')
 1320 FORMAT (1H ,'NEW EP=',E12.4,'    T=',D15.7)
 1330 FORMAT (1H ,'H IS TOO SMALL, H=',D12.5)
 1340 FORMAT (1H ,'WRONG CONTROL PARAMETERS')
C
      END
